How to work with AppEEARS Cloud Optimized GeoTIFF (COG) outputs¶

Summary

This tutorial demonstrates how to access AppEEARS Cloud Optimized GeoTIFF (COG) outputs. NASA's Application for Extracting and Exploring Analysis Ready Samples (AρρEEARS) has been migrated to NASA's Earthdata Cloud space located in AWS us-west 2. This enables the user working in the cloud instance deployed in AWS us-west 2 to access outputs direcly in the cloud using S3 link returned in the location header of the response. In this tutorial, we will walk through the process of submitting an area sample and accessing a Cloud Optimized GeoTIFF (COG) outputs from AppEEARS. The Dixie Fire, the second-largest fire in California history is used as an example in this tutorial. According to CalFire, the fire has started on July 13, 2021 and burned more than 963,276 acres. On August 18, the Dixie Fire merged with the Morgan Fire, which had been started by lightning August 12, close to Lassen National Park. The fire was one hundred percent contained by October 2021.

Requirements

  • Earthdata Login Authentication is required to access AppEEARS API and AppEEARS outpurs direcrly from an Amazon AWS bucket. See Requirements section in README.md.

Learning Objectives

  • Learn how to access AppEEARS Cloud Optimized GeoTIFF (COG) outputs

Tutorial Outline

  1. Setting Up
  2. Submit an area request in AppEEARS
  3. Extract the Direct S3 links
  4. Create a boto3 Refreshable Session
  5. Single COG File In-Region Direct S3 Access
  6. Multiple COG File In-Region Direct S3 Access
  7. Explore the EVI Time Series

1. Set up¶

Import the required packages.

In [1]:
import requests
import getpass, pprint, time, os, cgi, json
import geopandas 
import numpy
import datetime
import os
import io
import json
from urllib import parse
import requests
from netrc import netrc
from uuid import uuid4
import time
from pathlib import Path
from datetime import datetime, timezone
from botocore.client import Config
import rioxarray
import xarray
import hvplot.xarray
import holoviews
import geoviews
import rasterio 
from rasterio.plot import show
from rasterio.session import AWSSession
import s3fs
import pandas
import warnings
import sys
sys.path.append('../modules/')
import aws_session
warnings.filterwarnings('ignore')

In order to successfully run this tutorial, it is required to create a .netrc file in your home directory. The function _validate_netrc defined in aws_session checks if the netrc file with proper format exists in your home directory. If the netrc file does not exist, it will prompt to ask your Earthdata Login username and password and will create a netrc file. Please see the Prerequisites section in README.md.

In [2]:
# validate if netrc file is present else create one via the user / password prompt for the urs.earthdata.nasa.gov
aws_session._validate_netrc()

INFO : netrc file is setup for urs.earthdata.nasa.gov already ...

2. Submit an area request in AppEEARS¶

In this step, we are going to submit an area request with GeoTIFF as an output format. You can also submit this request using AppEEARS Graphic User Interface (GUI) and upload the JSON file provided in the repository (AppEEARS-Data-Resources/additional_file/Dixie-Fire-request.json). If you have your completed request, save your task_id to a variable, skip this step, and move to the next step of tutorial.

Assign the AρρEEARS API endpoint to a variable.

In [3]:
appeears_API_endpoint = 'https://appeears.earthdatacloud.nasa.gov/api/'

To access AppEEARS API a token is needed. This token is created using AppEEARS API endpoint, Earthdata Login credential stored in .netrc file, and requests.post function. This generated token will be added to the header. The header will be passed to requests.post function showing you are a qualifid user for submit and access a request.

In [4]:
urs = 'urs.earthdata.nasa.gov'
token_response = requests.post('{}login'.format(appeears_API_endpoint), auth = (netrc().authenticators(urs)[0],netrc().authenticators(urs)[2])).json() # Insert API URL, call login service, provide credentials & return json 
token = token_response['token']                      # Save login token to a variable
head = {'Authorization': 'Bearer {}'.format(token)}  # Create a header to store token information, needed to submit a request

Next, compile a JSON object with the request parameters. Dixie fire, started on July 13, 2021, but we extended the search query to two years to see the time series. The GeoJSON of Region of Interest(ROI) including Lassen National Park region, CA can be downloaded from the repository. For this tutorial, we are requesting _500m_16_days_EVI layer of MOD13A1.061 to see how Enhanced Vegetation Indices (EVI) varies before and after the fire event. Learn more about the MODIS Vegetation Indices 16-Day Version 6.1 product here. Below the AppEEARS search parameters are defined and the task JSON .

In [5]:
task_name = "Dixie Fire"
task_type = 'area'                  # Type of task, area or point
proj = 'geographic'                 # Set output projection 
outFormat = 'geotiff'               # Set output file format type
startDate = '01-01-2021'            # Start of the date range for which to extract data: MM-DD-YYYY
endDate = '12-31-2022'              # End of the date range for which to extract data: MM-DD-YYYY
ROI =  geopandas.read_file('../additional_files/DixieFire.geojson').to_json()
prodLayer = [{'layer': '_500m_16_days_EVI', 'product': 'MOD13A1.061'}]
In [6]:
task = {
    'task_type': task_type,
    'task_name': task_name,
    'params': {
         'dates': [
         {
             'startDate': startDate,
             'endDate': endDate
         }],
         'layers': prodLayer,
         'output': {
                 'format': {
                         'type': outFormat}, 
                         'projection': proj},
         'geo': json.loads(ROI),
    }
}

Next, submit the AppEEARS request using post function from requests library.

In [7]:
task_response = requests.post('{}task'.format(appeears_API_endpoint), json=task, headers=head).json()   # Post json to the API task service, return response as json
task_response                                                                  # Print task response
Out[7]:
{'task_id': 'a76dc9ab-867f-4cd0-bc26-a5eb7b54820e', 'status': 'pending'}

Save the task_id and wait until your request is processed and complete.

In [8]:
task_id = task_response['task_id']
task_id
Out[8]:
'a76dc9ab-867f-4cd0-bc26-a5eb7b54820e'
In [9]:
# Ping API until request is complete, then continue to Section 3
while requests.get('{}task/{}'.format(appeears_API_endpoint, task_id), headers=head).json()['status'] != 'done':
    print(requests.get('{}task/{}'.format(appeears_API_endpoint, task_id), headers=head).json()['status'])
    time.sleep(60)
print(requests.get('{}task/{}'.format(appeears_API_endpoint, task_id), headers=head).json()['status'])

queued
queued
processing
processing
done

3. Extract the Direct S3 links¶

Now that we have our outputs ready, we can get the bundle information for the files included in the outputs. If you submitted your request using AppEEARS GUI, assign your sample's task_id to the variable task_id below.

In [10]:
# task_id = 'a1714bef-bb86-4bd3-be54-bb84c01a11d8'

requests.get is used toget the bundle information. Below, bundle information for the first output file is printed. The bundle information includes s3_url in addition to the other information such as file_name, file_id, and file_type.
Each output file can be downloaded using the file_id and AppEEARS API endpoint. AppEEARS outputs are stored in an AWS bucket that can be accessed using S3_url.

In [11]:
bundle = requests.get('{}bundle/{}'.format(appeears_API_endpoint,task_id), headers=head).json()  # Call API and return bundle contents for the task_id as json

bundle['files'][0]
Out[11]:
{'sha256': 'ff02801f5d242355793a4942e705a07494d08c4671c19f6113f730e75d145fe3',
 'file_id': 'd44891bd-add8-415d-8092-0d7222e54b26',
 'file_name': 'MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2020353_aid0001.tif',
 'file_size': 341889,
 'file_type': 'tif',
 's3_url': 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2020353_aid0001.tif'}

Below, the S3 Links to Cloud Optimized GeoTIFF outputs and then S3 links for EVI layers are filted.

In [12]:
cog_urls = [f['s3_url'] for f in bundle['files'] if f['file_type'] == 'tif']
# cog_urls
In [13]:
EVI_cog_urls = [link for link in cog_urls if '500m_16_days_EVI' in link]
EVI_cog_urls
Out[13]:
['s3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2020353_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2021001_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2021017_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2021033_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2021049_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2021065_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2021081_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2021097_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2021113_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2021129_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2021145_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2021161_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2021177_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2021193_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2021209_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2021225_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2021241_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2021257_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2021273_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2021289_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2021305_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2021321_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2021337_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2021353_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2022001_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2022017_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2022033_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2022049_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2022065_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2022081_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2022097_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2022113_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2022129_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2022145_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2022161_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2022177_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2022193_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2022209_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2022225_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2022241_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2022257_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2022273_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2022289_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2022305_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2022321_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2022337_aid0001.tif',
 's3://appeears-output/a76dc9ab-867f-4cd0-bc26-a5eb7b54820e/MOD13A1.061_2020352_to_2022365/MOD13A1.061__500m_16_days_EVI_doy2022353_aid0001.tif']

In the same way, get the links to quality layers.

In [14]:
qual_cog_urls = [link for link in cog_urls if '500m_16_days_VI_Quality' in link]

4. Create a boto3 Refreshable Session¶

AppEEARS outputs are freely accessible from a cloud instance in us-west-2 region. In order to access our output files, a Boto3 session is needed. The Boto session stores the required configurations for an easy integration between Python and AWS services. Below, get_boto3_refreshable_session stored in aws_session will access your Earthdata login credentidals store in .netrc file and generate S3 credential by making a call to AppEEARS S3 credential endpoint, and create a boto3 session. This session will be auto-renewed as needed to prevent timeouts errors related to S3 credentials.

In [15]:
region_name = 'us-west-2'
s3_creds_endpoint = f"{appeears_API_endpoint}/s3credentials"
In [16]:
# Boto3 Session required by the rioxarray package 
boto3_session = aws_session.get_boto3_refreshable_session(s3_creds_endpoint, region_name)

INFO : Getting new temporary credentials from AppEEARS API https://appeears.earthdatacloud.nasa.gov/api//s3credentials...

Below, GDAL and AWS configuration and any environmner options are passed to rasterio.env(). When the Python context manager is entered all the configuration options are set and when the context is exited, configuration options are removed.

In [17]:
rio_env = rasterio.Env(
    AWSSession(boto3_session),
    GDAL_DISABLE_READDIR_ON_OPEN='EMPTY_DIR',
    GDAL_HTTP_COOKIEFILE=os.path.expanduser('~/cookies.txt'),
    GDAL_HTTP_COOKIEJAR=os.path.expanduser('~/cookies.txt')
)
rio_env.__enter__()
Out[17]:
<rasterio.env.Env at 0x7fa552504910>

5. Single COG File In-Region Direct S3 Access¶

COG datasets can be read using open_rasterio function from rioxarray library. The coordinates are generated automatically from the file’s geoinformation. Below, the first EVI S3 link is read as a xarray.DataArray. The data array is scaled and nodata values are masked by setting mask_and_scale to True. Chunk sizes can be provided or can be set to 'auto' to make a sensible chunk sized according to each dimension.

In [18]:
EVI_obj = rioxarray.open_rasterio(EVI_cog_urls[0], chunks = 'auto'  , mask_and_scale = True).squeeze('band', drop=True)

EVI_obj
Out[18]:
<xarray.DataArray (y: 412, x: 399)>
dask.array<getitem, shape=(412, 399), dtype=float32, chunksize=(412, 399), chunktype=numpy.ndarray>
Coordinates:
  * x            (x) float64 -121.8 -121.8 -121.8 ... -120.1 -120.1 -120.1
  * y            (y) float64 41.21 41.21 41.2 41.2 ... 39.51 39.51 39.5 39.5
    spatial_ref  int64 0
Attributes:
    AREA_OR_POINT:  Area
    units:          NA
xarray.DataArray
  • y: 412
  • x: 399
  • dask.array<chunksize=(412, 399), meta=np.ndarray>
    Array Chunk
    Bytes 642.14 kiB 642.14 kiB
    Shape (412, 399) (412, 399)
    Dask graph 1 chunks in 3 graph layers
    Data type float32 numpy.ndarray
    		399 412
    • x
      (x)
      float64
      -121.8 -121.8 ... -120.1 -120.1
      array([-121.797917, -121.79375 , -121.789583, ..., -120.147917, -120.14375 ,
       -120.139583])
    • y
      (y)
      float64
      41.21 41.21 41.2 ... 39.5 39.5
      array([41.210417, 41.20625 , 41.202083, ..., 39.50625 , 39.502083, 39.497917])
    • spatial_ref
      ()
      int64
      0
      crs_wkt :
      GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
      semi_major_axis :
      6378137.0
      semi_minor_axis :
      6356752.314245179
      inverse_flattening :
      298.257223563
      reference_ellipsoid_name :
      WGS 84
      longitude_of_prime_meridian :
      0.0
      prime_meridian_name :
      Greenwich
      geographic_crs_name :
      WGS 84
      grid_mapping_name :
      latitude_longitude
      spatial_ref :
      GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
      GeoTransform :
      -121.79999998908852 0.004166666666293395 0.0 41.21249999630797 0.0 -0.004166666666293395
      array(0)
    • x
      PandasIndex
      PandasIndex(Float64Index([-121.79791665575537, -121.79374998908908, -121.78958332242279,
              -121.78541665575649,  -121.7812499890902,  -121.7770833224239,
              -121.77291665575761, -121.76874998909132, -121.76458332242503,
              -121.76041665575873,
              ...
              -120.17708332256724, -120.17291665590095, -120.16874998923466,
              -120.16458332256836, -120.16041665590207, -120.15624998923577,
              -120.15208332256948, -120.14791665590319,  -120.1437499892369,
               -120.1395833225706],
             dtype='float64', name='x', length=399))
    • y
      PandasIndex
      PandasIndex(Float64Index([ 41.21041666297482,  41.20624999630853, 41.202083329642235,
               41.19791666297594,  41.19374999630965, 41.189583329643355,
               41.18541666297706,  41.18124999631077, 41.177083329644475,
               41.17291666297818,
              ...
               39.53541666312488, 39.531249996458584,  39.52708332979229,
                 39.522916663126, 39.518749996459704,  39.51458332979341,
               39.51041666312712, 39.506249996460824,  39.50208332979453,
               39.49791666312824],
             dtype='float64', name='y', length=412))
  • AREA_OR_POINT :
    Area
    units :
    NA

If the chunk shape and array are the same, the chunk selected for the dataset is bigger in size. you can set the chunk size manually.

In [19]:
EVI_obj.chunk(chunks={"x": 5, "y": 5})
Out[19]:
<xarray.DataArray (y: 412, x: 399)>
dask.array<rechunk-merge, shape=(412, 399), dtype=float32, chunksize=(5, 5), chunktype=numpy.ndarray>
Coordinates:
  * x            (x) float64 -121.8 -121.8 -121.8 ... -120.1 -120.1 -120.1
  * y            (y) float64 41.21 41.21 41.2 41.2 ... 39.51 39.51 39.5 39.5
    spatial_ref  int64 0
Attributes:
    AREA_OR_POINT:  Area
    units:          NA
xarray.DataArray
  • y: 412
  • x: 399
  • dask.array<chunksize=(5, 5), meta=np.ndarray>
    Array Chunk
    Bytes 642.14 kiB 100 B
    Shape (412, 399) (5, 5)
    Dask graph 6640 chunks in 4 graph layers
    Data type float32 numpy.ndarray
    		399 412
    • x
      (x)
      float64
      -121.8 -121.8 ... -120.1 -120.1
      array([-121.797917, -121.79375 , -121.789583, ..., -120.147917, -120.14375 ,
       -120.139583])
    • y
      (y)
      float64
      41.21 41.21 41.2 ... 39.5 39.5
      array([41.210417, 41.20625 , 41.202083, ..., 39.50625 , 39.502083, 39.497917])
    • spatial_ref
      ()
      int64
      0
      crs_wkt :
      GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
      semi_major_axis :
      6378137.0
      semi_minor_axis :
      6356752.314245179
      inverse_flattening :
      298.257223563
      reference_ellipsoid_name :
      WGS 84
      longitude_of_prime_meridian :
      0.0
      prime_meridian_name :
      Greenwich
      geographic_crs_name :
      WGS 84
      grid_mapping_name :
      latitude_longitude
      spatial_ref :
      GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
      GeoTransform :
      -121.79999998908852 0.004166666666293395 0.0 41.21249999630797 0.0 -0.004166666666293395
      array(0)
    • x
      PandasIndex
      PandasIndex(Float64Index([-121.79791665575537, -121.79374998908908, -121.78958332242279,
              -121.78541665575649,  -121.7812499890902,  -121.7770833224239,
              -121.77291665575761, -121.76874998909132, -121.76458332242503,
              -121.76041665575873,
              ...
              -120.17708332256724, -120.17291665590095, -120.16874998923466,
              -120.16458332256836, -120.16041665590207, -120.15624998923577,
              -120.15208332256948, -120.14791665590319,  -120.1437499892369,
               -120.1395833225706],
             dtype='float64', name='x', length=399))
    • y
      PandasIndex
      PandasIndex(Float64Index([ 41.21041666297482,  41.20624999630853, 41.202083329642235,
               41.19791666297594,  41.19374999630965, 41.189583329643355,
               41.18541666297706,  41.18124999631077, 41.177083329644475,
               41.17291666297818,
              ...
               39.53541666312488, 39.531249996458584,  39.52708332979229,
                 39.522916663126, 39.518749996459704,  39.51458332979341,
               39.51041666312712, 39.506249996460824,  39.50208332979453,
               39.49791666312824],
             dtype='float64', name='y', length=412))
  • AREA_OR_POINT :
    Area
    units :
    NA

To view the values in the data array, you can use load function.

In [20]:
EVI_obj.load()
Out[20]:
<xarray.DataArray (y: 412, x: 399)>
array([[ 0.4951    ,  0.5316    ,  0.5316    , ...,  0.0829    ,
         0.0829    ,  0.0339    ],
       [ 0.4951    ,  0.4951    ,  0.5316    , ...,  0.1221    ,
         0.1343    ,  0.1433    ],
       [ 0.4894    ,  0.3274    ,  0.3274    , ...,  0.1858    ,
         0.1646    ,  0.1646    ],
       ...,
       [        nan, -0.0211    , -0.0211    , ...,  0.2356    ,
         0.1786    ,  0.0675    ],
       [ 0.0751    , -0.0254    , -0.0211    , ...,  0.1785    ,
         0.2066    ,  0.1476    ],
       [ 0.017     , -0.0402    ,         nan, ...,  0.1785    ,
         0.1785    ,  0.21519999]], dtype=float32)
Coordinates:
  * x            (x) float64 -121.8 -121.8 -121.8 ... -120.1 -120.1 -120.1
  * y            (y) float64 41.21 41.21 41.2 41.2 ... 39.51 39.51 39.5 39.5
    spatial_ref  int64 0
Attributes:
    AREA_OR_POINT:  Area
    units:          NA
xarray.DataArray
  • y: 412
  • x: 399
  • 0.4951 0.5316 0.5316 0.5316 0.5316 ... 0.1971 0.1785 0.1785 0.2152
    array([[ 0.4951    ,  0.5316    ,  0.5316    , ...,  0.0829    ,
         0.0829    ,  0.0339    ],
       [ 0.4951    ,  0.4951    ,  0.5316    , ...,  0.1221    ,
         0.1343    ,  0.1433    ],
       [ 0.4894    ,  0.3274    ,  0.3274    , ...,  0.1858    ,
         0.1646    ,  0.1646    ],
       ...,
       [        nan, -0.0211    , -0.0211    , ...,  0.2356    ,
         0.1786    ,  0.0675    ],
       [ 0.0751    , -0.0254    , -0.0211    , ...,  0.1785    ,
         0.2066    ,  0.1476    ],
       [ 0.017     , -0.0402    ,         nan, ...,  0.1785    ,
         0.1785    ,  0.21519999]], dtype=float32)
    • x
      (x)
      float64
      -121.8 -121.8 ... -120.1 -120.1
      array([-121.797917, -121.79375 , -121.789583, ..., -120.147917, -120.14375 ,
       -120.139583])
    • y
      (y)
      float64
      41.21 41.21 41.2 ... 39.5 39.5
      array([41.210417, 41.20625 , 41.202083, ..., 39.50625 , 39.502083, 39.497917])
    • spatial_ref
      ()
      int64
      0
      crs_wkt :
      GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
      semi_major_axis :
      6378137.0
      semi_minor_axis :
      6356752.314245179
      inverse_flattening :
      298.257223563
      reference_ellipsoid_name :
      WGS 84
      longitude_of_prime_meridian :
      0.0
      prime_meridian_name :
      Greenwich
      geographic_crs_name :
      WGS 84
      grid_mapping_name :
      latitude_longitude
      spatial_ref :
      GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
      GeoTransform :
      -121.79999998908852 0.004166666666293395 0.0 41.21249999630797 0.0 -0.004166666666293395
      array(0)
    • x
      PandasIndex
      PandasIndex(Float64Index([-121.79791665575537, -121.79374998908908, -121.78958332242279,
              -121.78541665575649,  -121.7812499890902,  -121.7770833224239,
              -121.77291665575761, -121.76874998909132, -121.76458332242503,
              -121.76041665575873,
              ...
              -120.17708332256724, -120.17291665590095, -120.16874998923466,
              -120.16458332256836, -120.16041665590207, -120.15624998923577,
              -120.15208332256948, -120.14791665590319,  -120.1437499892369,
               -120.1395833225706],
             dtype='float64', name='x', length=399))
    • y
      PandasIndex
      PandasIndex(Float64Index([ 41.21041666297482,  41.20624999630853, 41.202083329642235,
               41.19791666297594,  41.19374999630965, 41.189583329643355,
               41.18541666297706,  41.18124999631077, 41.177083329644475,
               41.17291666297818,
              ...
               39.53541666312488, 39.531249996458584,  39.52708332979229,
                 39.522916663126, 39.518749996459704,  39.51458332979341,
               39.51041666312712, 39.506249996460824,  39.50208332979453,
               39.49791666312824],
             dtype='float64', name='y', length=412))
  • AREA_OR_POINT :
    Area
    units :
    NA

Now the data are loaded, lets quickly visualize the first EVI file using hvplot.image function.

In [21]:
EVI_obj.hvplot.image(x='x', y='y', rasterize=True, title= EVI_cog_urls[0].split('/')[-1], colorbar=True, cmap='YlGnBu').opts(clim=(0, 0.8))
Out[21]:

6. Multiple COG File In-Region Direct S3 Access¶

Next, lets access the EVI time series. Here, the EVI time series is being read and concadenated to a single xarray.DataArray using the concat function.

To do this, we first build a list of times from the S3 URLs and convert that to an xarray variable.

In [22]:
time_list = []
for obs in range(len(EVI_cog_urls)):
    doy = EVI_cog_urls[obs].split('_doy')[1].split('_aid')[0]
    doy_date = datetime.strptime(doy, '%Y%j').strftime('%d-%m-%Y')
    time_list.append(doy_date)
time = xarray.Variable('time', time_list)
In [23]:
chunks=dict(band=1, x=100, y=100)
EVI_series = xarray.concat([rioxarray.open_rasterio(f, chunks = 'auto', mask_and_scale = True).squeeze('band', drop=True) for f in EVI_cog_urls], dim=time)
EVI_series = EVI_series.rename('EVI')
EVI_series
Out[23]:
<xarray.DataArray 'EVI' (time: 47, y: 412, x: 399)>
dask.array<concatenate, shape=(47, 412, 399), dtype=float32, chunksize=(1, 412, 399), chunktype=numpy.ndarray>
Coordinates:
  * x            (x) float64 -121.8 -121.8 -121.8 ... -120.1 -120.1 -120.1
  * y            (y) float64 41.21 41.21 41.2 41.2 ... 39.51 39.51 39.5 39.5
    spatial_ref  int64 0
  * time         (time) <U10 '18-12-2020' '01-01-2021' ... '19-12-2022'
Attributes:
    AREA_OR_POINT:  Area
    units:          NA
xarray.DataArray
'EVI'
  • time: 47
  • y: 412
  • x: 399
  • dask.array<chunksize=(1, 412, 399), meta=np.ndarray>
    Array Chunk
    Bytes 29.47 MiB 642.14 kiB
    Shape (47, 412, 399) (1, 412, 399)
    Dask graph 47 chunks in 189 graph layers
    Data type float32 numpy.ndarray
    		399 412 47
    • x
      (x)
      float64
      -121.8 -121.8 ... -120.1 -120.1
      array([-121.797917, -121.79375 , -121.789583, ..., -120.147917, -120.14375 ,
       -120.139583])
    • y
      (y)
      float64
      41.21 41.21 41.2 ... 39.5 39.5
      array([41.210417, 41.20625 , 41.202083, ..., 39.50625 , 39.502083, 39.497917])
    • spatial_ref
      ()
      int64
      0
      crs_wkt :
      GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
      semi_major_axis :
      6378137.0
      semi_minor_axis :
      6356752.314245179
      inverse_flattening :
      298.257223563
      reference_ellipsoid_name :
      WGS 84
      longitude_of_prime_meridian :
      0.0
      prime_meridian_name :
      Greenwich
      geographic_crs_name :
      WGS 84
      grid_mapping_name :
      latitude_longitude
      spatial_ref :
      GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
      GeoTransform :
      -121.79999998908852 0.004166666666293395 0.0 41.21249999630797 0.0 -0.004166666666293395
      array(0)
    • time
      (time)
      <U10
      '18-12-2020' ... '19-12-2022'
      array(['18-12-2020', '01-01-2021', '17-01-2021', '02-02-2021', '18-02-2021',
       '06-03-2021', '22-03-2021', '07-04-2021', '23-04-2021', '09-05-2021',
       '25-05-2021', '10-06-2021', '26-06-2021', '12-07-2021', '28-07-2021',
       '13-08-2021', '29-08-2021', '14-09-2021', '30-09-2021', '16-10-2021',
       '01-11-2021', '17-11-2021', '03-12-2021', '19-12-2021', '01-01-2022',
       '17-01-2022', '02-02-2022', '18-02-2022', '06-03-2022', '22-03-2022',
       '07-04-2022', '23-04-2022', '09-05-2022', '25-05-2022', '10-06-2022',
       '26-06-2022', '12-07-2022', '28-07-2022', '13-08-2022', '29-08-2022',
       '14-09-2022', '30-09-2022', '16-10-2022', '01-11-2022', '17-11-2022',
       '03-12-2022', '19-12-2022'], dtype='<U10')
    • x
      PandasIndex
      PandasIndex(Float64Index([-121.79791665575537, -121.79374998908908, -121.78958332242279,
              -121.78541665575649,  -121.7812499890902,  -121.7770833224239,
              -121.77291665575761, -121.76874998909132, -121.76458332242503,
              -121.76041665575873,
              ...
              -120.17708332256724, -120.17291665590095, -120.16874998923466,
              -120.16458332256836, -120.16041665590207, -120.15624998923577,
              -120.15208332256948, -120.14791665590319,  -120.1437499892369,
               -120.1395833225706],
             dtype='float64', name='x', length=399))
    • y
      PandasIndex
      PandasIndex(Float64Index([ 41.21041666297482,  41.20624999630853, 41.202083329642235,
               41.19791666297594,  41.19374999630965, 41.189583329643355,
               41.18541666297706,  41.18124999631077, 41.177083329644475,
               41.17291666297818,
              ...
               39.53541666312488, 39.531249996458584,  39.52708332979229,
                 39.522916663126, 39.518749996459704,  39.51458332979341,
               39.51041666312712, 39.506249996460824,  39.50208332979453,
               39.49791666312824],
             dtype='float64', name='y', length=412))
    • time
      PandasIndex
      PandasIndex(Index(['18-12-2020', '01-01-2021', '17-01-2021', '02-02-2021', '18-02-2021',
       '06-03-2021', '22-03-2021', '07-04-2021', '23-04-2021', '09-05-2021',
       '25-05-2021', '10-06-2021', '26-06-2021', '12-07-2021', '28-07-2021',
       '13-08-2021', '29-08-2021', '14-09-2021', '30-09-2021', '16-10-2021',
       '01-11-2021', '17-11-2021', '03-12-2021', '19-12-2021', '01-01-2022',
       '17-01-2022', '02-02-2022', '18-02-2022', '06-03-2022', '22-03-2022',
       '07-04-2022', '23-04-2022', '09-05-2022', '25-05-2022', '10-06-2022',
       '26-06-2022', '12-07-2022', '28-07-2022', '13-08-2022', '29-08-2022',
       '14-09-2022', '30-09-2022', '16-10-2022', '01-11-2022', '17-11-2022',
       '03-12-2022', '19-12-2022'],
      dtype='object', name='time'))
  • AREA_OR_POINT :
    Area
    units :
    NA

7. Explore the EVI Time Series¶

Below, a dynamic plot is created to visually look at the variation in EVI through out these two years. You can manually select the dates from the dropdown next to the map. The latest observation before the fire event is on July 12, 2021. and the next observation is on July 28, 2021 which is the first observation after the fire event. A sudden decrease in EVI after July 12, 2021 is visually noticable and the area that vegetation is being removed grows spatially as you look through next observations. If you look at the observattions of the next year, you still can identify the scar left by fire.

In [24]:
map = EVI_series.hvplot.image(x='x', y='y',groupby='time', rasterize=True, colorbar=True, frame_width= 800, frame_height= 500, cmap='YlGnBu').opts(clim=(0, 0.8))
map
Out[24]:

You can also create a plot showing the EVI time series for a specific latittude and longitude.

In [25]:
EVI_series.sel(x=-121.260, y=40.330,method='nearest').hvplot.line(x='time', y='EVI', rot=90, frame_width= 1000, frame_height= 300, fontscale=1.5)
Out[25]:

Now, combine the time series map and the plot. Below,a dynamic map is created that shows the EVI time series as you move your mouse to a different location on the map. You can update the base map by selecting the date from the dropdown menu on the far right side. This is an easy way to explore your data before further processing.

In [26]:
# Stream of X and Y positional data
posxy = holoviews.streams.PointerXY(source=map, x=-121.4, y=40.03) 

# Function to build a time series using the mouse hover positional information retrieved from the map 
def point_spectra(x,y):
    return EVI_series.sel(x=x,y=y,method='nearest').hvplot.line(x='time', y='EVI', rot=90, frame_width= 800, frame_height= 500, fontscale=1.5) 

# Define the Dynamic Maps
point_dmap = holoviews.DynamicMap(point_spectra, streams=[posxy])

# Plot the Map and Dynamic Map side by side
map + point_dmap
Out[26]:

Finally, exit the context manager to remove the configuration options.

In [27]:
# Exit our context
rio_env.__exit__()

Contact Info:¶

Email: LPDAAC@usgs.gov
Voice: +1-866-573-3222
Organization: Land Processes Distributed Active Archive Center (LP DAAC)¹
Website: https://lpdaac.usgs.gov/
Date last modified: 02-16-2022

¹Work performed under USGS contract G15PD00467 for NASA contract NNG14HH33I.